Revisiting algorithms for generating surrogate time series 
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The method of surrogates is one of the key concepts of nonhnear data analysis. Here, we demon- 
strate that commonly used algorithms for generating surrogates often fail to generate truly linear 
time series. Rather, they create surrogate realizations with Fourier phase correlations leading to 
non-detections of nonlinearities. We argue that reliable surrogates can only be generated, if one 
tests separately for static and dynamic nonlinearities. 
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Introduction. The method of surrogates [T] is one of 
the the key concepts of nonlinear data analysis, which 
allows to test for weak nonlinearities in data sets in a 
model-independent way. The basic idea of this approach 
is to compute statistics sensitive to nonlinearities for the 
original data set and for an ensemble of so-called surro- 
gate data sets, which mimic the linear properties of the 
original data. If the computed measure for the original 
data is significantly different from the values obtained for 
the set of surrogates, one can infer that the data contain 
nonlinearities. 

Linearity means in this case that all the structure in the 
time series is contained in the autocorrelation function, 
or equivalently, in the Fourier power spectrum. Thus the 
time series yt can be modeled e.g. by an autoregressive 
(AR) model described by j/t = oq + YX=i O'kUt-k + cret 
(a: coefficients, e: white noise) or a more general ARMA 
models also including moving averages (MA). Nonlinear- 
ity thus refers to all those structures in data sets that are 
not captured by the power spectrum. 
Since its introduction the method of surrogates has found 
numerous applications in many fields of research ranging 
from geophysical and physiological time series analysis 
[2] to econophysics (3j, astrophysics and cosmology 

The most commonly used methods for generating sur- 
rogates include Fourier-transformed (FT), amplitude ad- 
justed Fourier-transformed (AAFT) [T] and iterative am- 
plitude adjusted Fourier-transformed (lAAFT) surro- 
gates [7]- With FT surrogates, which are generated by 
using the Fourier amplitudes of the original time series 
and a set of random phases, one tests whether the time 
series is compatible with a linear Gaussian process. In 
other words, any significant deviation of the original data 
from its surrogates indicates the presence of higher order 
temporal correlations, i.e. the presence of dynamic non- 
linearities. 

The null hypothesis can be generalized to cases where 
the data set is non-Gaussian. Here, the hypothesis is 
that the original data are from a nonlinear or linear 
stochastic process that has undergone a static nonlinear 



transformation. The AAFT algorithm was developed to 
generate this kind of surrogates. The original time series 
is rendered Gaussian by a rank-ordered remapping of the 
values onto a Gaussian distribution. For this Gaussian 
time series which follows the measured time evolution, a 
FT surrogate time series is generated. The final surro- 
gate data are obtained by rank-ordered remapping of the 
FT surrogate onto the distribution of the original data. 
It was shown that the final remapping step in the AAFT 
algorithm can lead to a whitening of the power spectrum 
7J . This shortcoming was overcome by the lAAFT algo- 
rithm, which consists of the following iteration scheme. 
(1) One starts with a random shuffle {?/^} of the orig- 
inal data {2/n}, which is Fourier transformed yielding 
the Fourier amplitudes |>S'*(fc)| and Fourier phases 0*(A:) 
. (2) The Fourier amplitudes [^'(fc)! of the time series 
{vli} are then replaced by those for the original data 
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Fourier (j)^{k) components are kept. An inverse Fourier 
transformation leads to a time series with desired power 
spectrum {z^}. (3) To make the surrogate time series 
have the same amplitude distribution {?/„} in real space 
as the original data, {2^} is remapped in a rank or- 
dered way onto {y^j} leading to {yj,^^}. As in AAFT 
the final remapping step changes the power spectrum, so 
that steps (2) and (3) have to be repeated until conver- 
gence to the correct power spectrum is achieved. AAFT 
and lAAFT strive for a proper reproduction of the lin- 
ear properties and the amplitude distribution of the time 
series. All three algorithms have in common that they in- 
tend to translate the absence of nonlinear temporal corre- 
lations into the constraint of uncorrelated and uniformly 
distributed phases in the Fourier representation of the 
data. While the randomness of the phases remains (by 
construction) exactly preserved for FT surrogates, it has 
never been studied, how the remapping step (in AAFT) 
and the iteration steps (in lAAFT) affect the randomness 
of the phases in the resulting surrogate time series. 
In this Letter we investigate in detail how well this much 
more important constraint of the absence of nonlinear 
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correlations in the surrogates is fulfilled for AAFT and 
lAAFT surrogates and study consequences for the out- 
come of tests on nonlinearity for observational data. 

Data Set and Preprocessing. All following investiga- 
tions are based on two - quite distinct - data sets. The 
first one is an X-ray observation of the narrow-line Seyfert 
1 galaxy Mrk 766 taken with the XMM Newton satellite 
[8]. For light curves from active galactic nuclei (AGN) 
one investigates the linear and nonlinear temporal cor- 
relations to infer more information about the physical 
processes in the innermost regions of these compact ac- 
creting objects. The time series originates from observa- 
tions of Mrk 766 in 2005. We selected the revolution 999 
from May, 23rd, 2005 that took 95 ks. We used the data 
retrieved from the XMM public archive and only relied 
on data from the PN camera [9] due to its superior sta- 
tistical quality. Source counts were accumulated from a 
rectangular box of 27 x 26 RAW pixels (1 RAW pixel ^ 
4.1 arcsec) around the position of the source. Background 
data were extracted from a similar, source free region on 
the chip. After rejection of times affected by high back- 
ground a light curve was extracted in the 0.3—10 keV 
energy band, background subtracted and binned to 50 s 
time resolution. 

The second data set are the day-to-day returns of the 
Dow Jones (D J) industrial index covering the period from 
26th of May 1896 until 8th of June 2012. From the pub- 
licly available (www.djaverages.com) closing prices p{t) 
of each trading day the day-to-day returns r{t) are de- 
termined by considering the change of the logarithm of 
the stock price rAt{t) = ^ lii(p(t — At)), where 

At — 1 day. Both data sets (see Fig. IT]) arc particularly 
well suited to be used for testing concepts of nonlinear 
time series analysis, because they represent scalar, real 
time series from a sufficiently complex system, where the 
mere detection of nonlinearities already allows to discrim- 
inate between different classes of models. In the case of 
the the AGN-data comparably simple intrinsically linear 
models for light curve variations like e.g. global disk os- 
cillations models [TD] can be safely ruled out. The detec- 
tion of signatures of nonlinear determinism in financial 
data shows that the stock market is neither completely 
governed by stochastic processes nor can it be purely de- 
scribed by linear processes plus additive noise. 

Analyses and Results. For both time series we gener- 
ate 200 reahzations of FT, AAFT and lAAFT surrogates 
each. Note that in the FT case the random phase hypoth- 
esis for the surrogates is doubtlessly fulfilled. To assess 
the presence or absence of nonlinear correlations in the 
surrogates we make use of the representation of Fourier 
phases in so-called phase maps [Ilj. A phase map is de- 
fined as a two-dimensional set of points G = {4>k, 0/c+a} 
where is the phase of the fc*'' mode of the Fourier 
transform and A a phase delay. If the phases were taken 
from a random, uniform and uncorrelated distribution, 
the phase maps G would be a random two-dimensional 
distribution of points in the interval [— 7r;7r]. Any sig- 
nificant deviation from this random distribution already 
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FIG. 1: Light curve of the observation of Mrk 766 in 50 s 
time resolution (above) and the day-to-day returns of the Dow 
Jones Industrial index (below). The time series consists of 
1540 and 29070 points respectively. 



points towards the presence of nonlinearities in the data. 
In Fig. |2]we show the phase maps for the AGN case for 
a delay A = 1 for one AAFT and one lAAFT realiza- 
tion. For the Dow Jones data we show the respective 
plots for A = 3. IFor both types of times series and for 
both classes of surrogates the phase maps are far from 
being random. To quantify the degree of correlation be- 
tween the phases <j) and A we calculate the correlation 
coefficient c(A) for delays up to A = 10. In Fig. |3]we 
show c(A) for the three classes of surrogates and both 
time series under study. The jitter around zero of the 
FT surrogates (red) represents - by construction - the 
statistical fiuctuations of c(A) for a limited number of 
phases. For both examples, however, these fluctuations 
around zero are significantly larger for the AAFT surro- 
gates (blue) and lAAFT surrogates (black) for a number 
of values of A. This results in up to 138 (I A AFT) and 
110 (AAFT) surrogate realizations leading to c- values ly- 
ing outside the 3 a confidence region around zero. Since 
the only difference between FT and AAFT is the remap- 
ping step after the inverse Fourier transformation, we 
can immediately infer that this step is responsible for in- 
ducing the phase correlations. The deviations from zero 
for the correlation coefficient become in both cases more 
pronounced for the lAAFT surrogates. For the AGN 
light curve one can even observe that for A = 1 the re- 
gion around zero for c(A) is thinned out meaning that 
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it is rather a rule than an exception that |c(A)| becomes 
larger than the normal statistical fluctuations. The pres- 
ence of phase correlations in the majority of surrogate 
realizations proves already in a strict mathematical sense 
that those realizations are in fact nonlinear. 
To address the relationship between phase correlations 
and measures for nonlinearity we calculate the nonlin- 
ear prediction error (NPLE) as described in [H]. The 
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FIG. 2: Phase maps for A = 1 for one AAFT (left) and one 
lAAFT (right) surrogate realization for the AGN time series 
(above). Below the phase maps for A = 3 for one AAFT 
(left) and one I A AFT (right) surrogate realization for the 
Dow Jones data are shown. 

NLPE has been shown to be a robust measure with a 
good overall performance [13]. The calculation of the 
NLPE relies on the representation of the time series in 
an artificial phase space, which is obtained using the 
method of delay coordinates [M]. This is accomplished 
by using time delayed versions of the observed time series 
as coordinates for the embedding space. The multivari- 
ate vectors in the d-dimensional space are expressed by 
Vn = {yn-{d-i)T, yn-(d-2)r, ■■■,yn) , whcrc T is the delay 
time and j/„ denotes the value of the (discretized) time 
series at time step n. The comparison of the predicted 
behavior of the time series based on the local neighbors 
with the real trajectory of the system leads to the defi- 
nition of the NLPE -0 as 



= ^Pid,T,T,N) 
1 



(1) 



(A/-r-(d-l)r)^ 



M-l-T 
ri=(d-l)T 



where F is a locally constant predictor, M is the length 
of the time series, and T is the lead time. The predic- 



FIG. 3: Cross-correlation between 0(fc) and (j){k + A) for A = 
1, . . . , 10 for 100 reahzations of FT (red), AAFT (blue) and 
lAAFT (black) surrogates for the Mrk 766 (above) and the 
Dow Jones (below) time series. The black (red) crosses denote 
the values for c(A) for the original time series (remapped 
original time series) for comparison. 



tor F is calculated by averaging over future values of the 
N {N = d+1) nearest neighbors in the delay coordinate 
representation. 

We found that for T > b %p remains rather constant, thus 
a value of T = 5 was used for this study. The dimension 
of the embedding space d and the delay time r have to 
be set appropriately. Since the AGN time series of this 
study consist only of 1539 data points , we use a low 
embedding dimension c? = 3, so that the vectors y„ are 
not too sparsely distributed in the embedding space. To 
allow for direct comparison we also used d = 3 for the 
much longer time series of the Dow Jones. 
The delay time r is determined by using the criterion of 
zero crossing of the autocorrelation function or the first 
minimum of the autocorrelation function or mutual in- 
formation [IS], which lies around r = 10 ks for Mrk 766 
and T = 2 days respectively. To also investigate the ro- 
bustness of our results with respect to variations of r, 
we calculate the NPLE for a set of four different delay- 
times, r = 7.5,8.75,10.0 and 11.25 ks (AGN data) and 
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T = 2,3,4: and 5 days (DJ data) and analyze the results 
as a function of r. 

First, we investigated the relation between the phase cor- 
relations c(A) and the NLPE iIj{t) for both time series. 
While no obvious trends between linear correlations in 
the phase maps and the nonlinear prediction error could 
be identified for the DJ data, we found remarkable anti- 
correlations for the AAFT (-0.68) and lAAFT (-0.94) 
between c(A — 1) and ^/'(t — lOfcs) for the Mrk 766 
data. The latter value means a nearly perfect correlation 
indicating that the NLPE calculated for lAAFT surro- 
gates is mostly determined by Fourier phase correlations 
for A = 1. We note here that this kind of analysis has 
the potential to shed more light on the effects of phase 
correlations on the shape of the point distribution in em- 
bedding space to ultimately get more insight into the 
meaning of Fourier phases for nonlinear time series. For 
this study, however, it is already sufficient to show that 
there are significant phase correlations which affect the 
calculation of nonlinear statistics. 

The results for a standard test on nonlinearities using the 
NLPE as test statistics are summarized in Fig. |4j We 
plot the significances S{ip), 

^^^-j '^original surrogate 

as a function of r. The mean {ijj) surrogate and the stan- 
dard deviation ct^^^^^^^^j^ are derived by summing over 
N = 200 realizations of the respective set of surro- 
gates. Obviously, the propagation of the phase corre- 
lations found for the AAFT and lAAFT surrogates to 
the calculation of the NLPE eventually leads to a non- 
detection of nonlinearities. For both time series the sig- 
nificances for the AAFT and lAAFT test are well below 
the 3 a detection criterion, whereas the test on dynamic 
nonlinearities using the remapped time series and FT 
surrogates yields a highly significant detection of nonlin- 
earities with S{ip) reaching eight for Mrk 766 and six for 
the Dow Jones . These results remain stable for a range 
of delay times. 

Conclusions. Using the examples of an AGN light 
curve and the Dow Jones day-to-day returns we demon- 
strated that both the AAFT and lAAFT algorithm can 
generate surrogate data sets with phase correlations. 
Thus, these time series have to be considered as non- 
linear. We further showed that those phase correlations 
propagate into the calculation of nonlinear statistics like 
the NLPE. Hence, nonlinearities in time series may re- 
main undetected due to the presence of them in a number 
of surrogate realizations against which the original data 
are compared. The wrong outcome of the surrogate test 
leads then, in turn, to a wrong modeling of the complex 
underlying system. In our case the AAFT and lAAFT 
surrogate test would suggest that the X-ray variability 
of Mrk 766 is compatible with a global disk oscillations 
model, which is clearly ruled out by the FT surrogate 
test. Similarly, the detection of weak dynamic nonlinear- 
ities (only with FT-surrogates) proves once again that 
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FIG. 4: Significances S{ip) for nonlinear temporal correlations 
for the set of light curves as determined with NLPE and as a 
function of the delay time r for Mrk 766 (left) and the Dow 
Jones (right). Red: results obtained using the remapped time 
series and FT surrogates. Black: results obtained with the 
original data and the ITAAFT algorithm. Blue: results for 
the original data and the AAFT algorithm. The dashed line 
indicates the 3 a detection limit. 



the returns are correlated with each other, which dis- 
proves one of the basic assumptions of the Black-Scholes 
model |16l I17j . where they are assumed to be random 
and independent from each other. Since the correlations 
are nonlinear, simple ARMA processes cannot be used as 
market models. Finally, the detected nonlinearities for 
the remapped data represent signatures that go beyond 
the so-called stylised facts (see e.g. [18]) of the fiuctua- 
tions of price indices, that mostly relate to the shape of 
the probability density derived from the original distri- 
bution of the return values, e. g. 'fat tails' etc. Market 
models must thus not only be able to reproduce these 
features but must also account for the intrinsic dynamic 
nonlinearities detected by means of FT-surrogate test. 
By analyzing a larger set of AGN light curves as well as 
other typical simulated nonlinear time series like e.g. the 
Lorenz system we convinced ourselves that the presence 
of phase correlations is a generic property of AAFT and 
lAAFT surrogates - rather independent of the underly- 
ing time series under study. 

Having outlined the consequences of surrogate tests on 
model building using the examples of an AGN light curve 
and the Dow Jones day -to-day returns it becomes ob- 
vious that it may be necessary that previous results ob- 
tained with AAFT and lAAFT have to be critically re- 
assessed. It may be that weak nonlinearities remained 
undetected when AAFT and lAAFT surrogates were 
used. As a consequence wrong physical models may have 
been found to be compatible with the data. 
In general, surrogate generating algorithms aiming at 
testing both static and dynamic nonlinearities cannot re- 
produce both the power spectrum and the amplitude dis- 
tribution while preserving the randomness of the Fourier 
phases so far. Thus, one has to test separately for static 
and dynamic nonlinearities. The latter can be achieved 
by using FT surrogates for remapped time series as a 
reliable and truly linear class of surrogates. Additional 
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comparably easy tests on the Gaussianity of the amph- 
tude distribution can exclude the presence of static non- 
linearities induced by a nonlinear rescaling. 
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